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Synthetic  Aperture  Sonar  Imaging  of  Simple 

Finite  Targets 

Steven  G.  Kargl,  Kevin  L.  Williams,  and  Eric  1.  Thorsos 

Abstract 

During  the  Sediment  Acoustics  Experiment  2004  (SAX04),  a  synthetic  aperture  sonar  (SAS)  was 
used  to  detect  simple  targets  that  were  either  proud  or  buried  below  a  water-sediment  interface,  where 
the  nominal  grazing  angle  of  incidence  from  the  SAS  to  the  point  above  a  buried  target  was  well  below 
the  critical  grazing  angle.  SAS  images  from  other  measurements  below  the  critical  angle  have  also 
produced  target  detections  of  buried  spheres  and  finite  cylinders.  Models  and  numerical  simulations  are 
developed  to  investigate  these  proud  and  buried  target  detections.  For  buried  targets,  the  simulations 
include  estimates  of  reverberation  from  the  rough  seafloor,  the  subcritical  penetration  through  the 
seafloor,  scattering  from  a  target,  and  propagation  back  to  the  SAS.  For  proud  targets,  the  simulations 
include  the  scattering  from  the  target  where  interaction  with  the  seafloor  is  included  through  simple 
ray  models.  The  simulations  used  environmental  and  material  parameters  measured  during  SAX04.  The 
environmental  measurements  include  profiles  of  small-scale  surface  roughness  and  superimposed  ripple 
structure.  The  SAS  simulations  and  model/measurement  comparisons  over  a  frequency  range  of  10-50 
kHz  further  support  scattering  from  sediment  ripple  structure  as  the  dominant  mechanism  for  subcritical 
penetration  in  this  range. 

Index  Terms 

Synthetic  aperture  sonar,  acoustic  scattering,  ripple,  signal-to-noise  ratio  (SNR) 

I.  Introduction 

Synthetic  aperture  sonar  (SAS)  systems  are  becoming  a  common  tool  in  underwater  acoustics, 
a  particular  example  being  systems  designed  to  detect  mines  in  shallow  water  [1],  [2].  In  such 
environments  the  water-sediment  interface  invariably  plays  a  role,  affecting  both  the  signal 
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from  the  target  being  detected  and  contributing  to  the  noise  that  masks  this  signal.  Thus, 
understanding  and  modeling  the  acoustic  processes  at  the  water- sediment  interface  and  simulating 
SAS  operation  in  the  presence  of  this  interface  are  important. 

The  primary  goals  of  this  work  are  to  summarize  the  modeling  techniques  used  to  quantitatively 
understand  SAS  imaging  of  simple  proud  and  buried  targets  and  compare  modeling  results  to  field 
data.  The  model  development  was  motivated  by  the  specific  problem  of  shallow  grazing  angle 
detection  of  buried  targets  at  frequencies  below  50  kHz.  However,  the  more  general  objective  of 
incorporating  the  environment  into  a  “single  realization”  modeling  capability  needed  to  generate 
SAS  images  was  equally  important.  Thus  several  problems  of  increasing  complexity  are  examined 
here  starting  with  SAS  imaging  of  freefield  targets  and  ending  with  the  buried  target  problem. 

The  field  data  used  for  comparisons  with  modeling  results  were  obtained  during  the  Sediment 
Acoustics  Experiment  in  2004  (SAX04)  [3].  SAX04  was  conducted  in  the  Gulf  of  Mexico  about 
1  km  offshore  near  Fort  Walton  Beach,  Florida  in  17-18  m  of  water.  The  SAS  measurements 
were  made  using  a  28  m  rail  system  [4]  deployed  on  the  sea  floor. 

Finite  element  modeling  (FEM)  can  also  be  applied  to  model  with  high  fidelity  scattering  from 
proud  and  buried  targets  (e.g.,  see  [5]).  However,  the  FEM  approach  is  computationally  intensive, 
and  therefore  its  application  to  SAS  imaging  that  involves  Fourier  synthesis  to  model  each 
transmission  would  be  especially  challenging.  The  approach  being  described  here  involves  a  more 
approximate  treatment  of  target  scattering,  but  has  the  advantage  of  much  better  computational 
efficiency  than  for  the  FEM  method. 

IF  Model  Descriptions 

Several  independent  models  have  been  combined  to  produce  a  single  model,  which  can  be 
used  to  generate  a  set  of  signals  for  SAS  processing.  This  section  describes  these  models  and 
some  of  the  underlying  assumptions  and  approximations. 

The  origin  of  the  coordinate  geometry  for  the  reverberation  and  penetration  models  lies  in 
the  mean  water- sediment  interface.  The  2;-axis  is  oriented  perpendicular  to  this  interface  such 
that  its  unit  vector,  2,  points  towards  the  air-water  interface.  The  x-axis  coincides  with  the  range 
coordinate  of  the  SAS.  The  y-axis  is  aligned  with  the  direction  of  travel  of  the  SAS,  sometimes 
called  the  along-track  direction  or  dimension.  Throughout  this  article,  the  y-axis  will  be  referred 
to  as  the  cross  range  direction  or  simply  as  cross  range.  For  a  given  signal,  the  locations  of 
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the  source,  receiver,  and  target  are  given  by  r^rc  =  {x^rcysrc,  Zsrc),  ^rcv  =  {xrcv,yrcv,  Zrcv),  and 
'f'tgt  =  {xtgt-,ytgt-,ztgt),  respectively.  For  the  surface  integrals  to  be  defined  shortly,  capital  bold 
symbols  denote  horizontal  vectors  (i.e.,  r  =  R-\-  zz). 

For  the  SAX04  experiment  and  the  simulations  to  be  discussed  below,  interaction  of  the 
acoustic  field  with  the  air-water  interface  does  not  occur  in  the  time  interval  of  interest.  Thus, 
the  model  assumes  semi-infinite  inviscid  water  above  a  semi-infinite  fluid  sediment.  The  sound 
speed  and  density  in  the  water  are  Ci  and  pi\  while  the  sound  speed,  density,  and  dimensionless 
loss  parameter  for  the  sediment  are  C2,  p2,  and  <52,  respectively  [6]. 


A.  Reverberation 


A  time  domain  reverberation  model  has  recently  been  developed,  based  on  first-order  pertur¬ 
bation  theory,  to  describe  sound  scattering  from  a  rough  seafloor.  A  derivation  of  this  model  is 
given  in  Appendix  I.  Several  approximations  are  employed  in  addition  to  first-order  perturbation 
theory:  the  stationary  phase  approximation  is  used  to  evaluate  the  zero-order  field  needed  for 
the  first-order  theory,  the  sound  speed  in  the  sediment  is  approximated  as  independent  of 
frequency,  the  attenuation  in  the  sediment  is  approximated  as  linearly  increasing  with  frequency, 
and  the  frequency  dependence  of  the  source  and  receive  beam  patterns  is  neglected.  The  latter 
approximation  is  reasonably  satisfied  for  our  measurements  because  the  horizontal  aperture  of 
the  source,  composed  of  four  transducers  in  a  compact  horizontal  array,  is  steadily  reduced  as 
the  frequency  is  increased  by  reducing  the  number  of  active  elements  from  four  down  to  one  in 
steps.  For  bistatic  backscattering  in  the  incident  vertical  plane  [see  the  discussion  near  (1.40)], 
this  model  reduces  to  the  following  surface  integral  over  the  mean  plane  of  the  water- sediment 
interface: 


Previ^'^rcv:  t') 


X 


^2 


(1) 


J  ^src| 

Here,  and  throughout  this  article,  the  assumed  time  convention  is  exp(— with  t  and  oj 
denoting  time  and  angular  frequency.  In  (1)  the  subscript  rev  denotes  quantities  related  to 
reverberation,  that  is,  sound  scattering  from  the  rough  water-sediment  interface.  Environmental 
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and  system  parameters  enter  into  the  formulation  through  the  eomplex  eoeffieient 

Arev{R)  =  f'oBsrcBrcv{[^  +  ril(0j)][l  +  rii(0s)] 

X  (1  -  n^p)  +  (p  -  1)[1  -  rn(0.)][l  -  Tnies)]  (2) 

X  sin  9i  sin  +  (1  -  1/ p)  [1  +  T n  (0*)]  [1  +  Fn  (6'J] 

X  cos^j  cos  9s}, 

where  Bgrc,  Brcv,  9i,  9s,  and  rii(6*)  depend  on  R.  The  source  and  receiver  response  functions 
evaluated  at  R  on  the  mean  surface  plane  are  Bsrc  and  Brcv,  respectively.  (The  beam  patterns  are 
the  square  of  the  response  functions.)  The  incident  and  scattered  grazing  angles  are  determined 
by  sin6*j  =  Zsrc/\R  —  ^srcl  and  sin6's  =  Zrcv/l^rcv  —  respectively.  The  plane-wave  reflection 
coefficient  at  grazing  angle  9  is 

psm9  -  {k^ -cos‘^9YB 

p  sm  6*  +  —  COS"*  9)  u ^ 

In  (3)  p  =  P2/P1  and  k  =  k2/ki  =  (1  + 1^2)1^,  where  v  =  C2/C1  is  the  real  index  of  refraction. 
(The  subscripts  on  rii(0)  denote  that  both  the  ingoing  and  outgoing  fields  are  in  the  water 
(medium  1),  in  contrast  to  the  transmission  coefficient  introduced  in  Sec.  II-B.)  Finally,  tq  =  1 
m  has  been  introduced  into  (2)  to  preserve  dimensional  consistency  in  (1).  In  (I),  h{R)  is  the  local 
height  of  the  water- sediment  interface  relative  to  the  mean  plane,  and  Psrc{t)  is  the  pressure  field 
at  1  m  from  the  source  on  the  transmit  beam  axis.  The  double  dot  notation  indicates  the  second 
derivative  with  respect  to  time.  Finally,  the  round-trip  time  delay  to  propagate  from  the  source  to 
the  point  R  on  the  mean  plane  and  then  back  to  the  receiver  is  Uev  =  {\rrcv  —  R\  +  \R—'rsrc\)/ci. 

B.  Penetration 

Thorsos  et  al.  [6]  proposed  rough  surface  scattering  as  a  mechanism  to  explain  the  significant 
penetration  of  an  incident  acoustic  field  into  sandy  sediments  below  the  critical  grazing  angle  9c  = 
cos“^(l/z/).  Subsequently,  Jackson  et  al.  [7]  demonstrated  that  ripples  are  the  primary  surface 
feature  needed  to  cause  this  increased  penetration.  Their  model  for  the  penetrating  acoustic 
pressure,  p,  is  expressed  in  a  baseband  representation,  i.e.,  the  signals  are  downshifted  by  the 
center  frequency.  In  the  present  case  the  signal  bandwidth  is  much  greater,  and  therefore  the 
baseband  representation  is  not  used,  but  otherwise  the  treatment  here  follows  closely  that  given 
in  [6].  Three  possible  contributions  are  contained  in  the  model:  a  zeroth-order  refracted  field. 
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PR^r,  t),  a  zeroth-order  evanescent  field,  Pe{r,  t),  and  a  first-order  surface  (i.e.,  interface)  scattered 
field,  P2s{r,  t).  Each  field  is  briefly  summarized  below  to  unify  the  notation  with  the  other  models 
presented  here  as  well  as  to  absorb  a  l/dvr  scaling  factor  in  (C15),  (C21),  and  (C25)  of  [6]  into  the 
amplitude  of  the  source  waveform.  This  results  in  the  same  normalization  for  the  reverberation 
and  penetration  models.  (In  [6]  the  source  waveform  normalization  was  taken  to  be  the  same 
as  the  point  source  Green’s  function  as  written  in  (Cl)  of  [6],  introducing  a  factor  of  l/An  in 
these  equations.) 

For  a  source  located  at  Vgrc  and  an  in-sediment  field  point  at  r/,  the  refracted  acoustic  field 
is  [6,  see  (C21)] 

PR{f'f-,i)  =  Bsrc,Rri2E{0iR)Psrc{t  —  tn) 

X  [\Rf-  Rsrc I  tan  92r  sin  9iR  (4) 

X  {\zsrc\/ sin^6'iR  -f  v\zf\/ siv?  92r)]~^^‘^ 

where  the  subscript  R  denotes  quantities  associated  with  the  refracted  path  connecting  the  source 
at  Tsrc  to  the  in-sediment  field  point  at  r/.  For  the  refracted  path,  9iR  denotes  the  in-water  grazing 
angle  at  the  interface,  and  6*2r  is  the  grazing  angle  of  the  refracted  ray  in  the  sediment  (medium 
2).  The  determination  of  9iR  will  be  described  below,  and  92r  is  related  to  9iR  by  Snell’s  law, 
cos  92r  =  E  cos  9iR,  where  the  loss  in  the  sediment  is  ignored  in  defining  92r-  The  use  of  absolute 
values  on  the  vertical  coordinates  ensures  that  (4)  is  invariant  to  the  2:-axis  convention.  Bsrc,R 
is  the  source  response  function  evaluated  at  the  refraction  point,  the  point  where  the  refracted 
path  intersects  the  mean  plane.  The  quantity  Ti2E{9iR)  is  a  real  transmission  coefficient  from 
medium  1  into  medium  2  defined  such  that  its  square  is  the  energy  transmission  coefficient.  It 
is  given  by 

ri2i<;(M  =  [i-|rii(MlT^',  (5) 

where  rii(6'ij?)  is  given  by  (3).  The  parameter  Ir  =  {diR  +  Kd2R)/ci  is  a  complex  time  delay 
where  the  imaginary  part  accounts  for  attenuation  in  the  sediment.  The  distance  from  the  source 
to  the  refraction  point  is  diR  =  \zsrc\/ and  d2R  =  |;2/|/ sin  6'2r  is  the  distance  from 
refraction  point  to  the  field  point  in  the  sediment. 

To  obtain  9iR,  we  begin  by  noting  that  the  horizontal  projections  of  the  distances  from 
the  source  to  the  refraction  point  and  from  the  refraction  point  to  the  field  point  are  given 
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by  k,,,!  cot  9iR  and  1 2:/ 1  cot  6^21?,  respectively.  These  must  combine  to  give  the  total  horizontal 
distance  between  the  source  and  field  point.  Thus,  we  have 

|-Rj  .^srcl  l^srcl  COt  OiJi  T  l^/l  ^2/?*  (b) 

It  readily  follows  from  (6)  that 

zjcot‘^62R  =  \Rf  -  Rsrcl"^  -  2\Rf  -  Rsrc\\Zsrc\  (V) 

X  cot  9iii  +  cot^  9iR. 


By  using  the  identity 

[l  +  (1  -  Z/^)  cot^  9iR]  cot^  92r  =  cot^  9iR,  (8) 

cot  6*2r  can  be  eliminated  from  (7)  and  the  following  quartic  equation  for  cot  9iR  is  obtained: 

4 

bj  cot^  9iR  =  0,  (9) 

j=0 

where 


^0  1"^/  -^srcl  •)  ^1  ^1"^/  -^src 1 1 ^src 1 5 

b2  =  + fi\Rf  -  R.rcl"^  -  u'^zj,  (10) 

63  2,fJ,\R,j^  1 1 I ,  (14  9'Zsrci 

with  /i  =  1  —  z/^.  Only  the  real,  positive  root  of  (9)  is  retained. 

For  the  evanescent  contribution,  we  define  the  incident  grazing  angle  9ie  to  be  the  grazing 
angle  at  the  point  on  the  interface  directly  above  the  field  point  in  the  sediment.  (The  subscript 
e  denotes  quantities  associated  with  the  evanescent  field.)  Thus, 

9ie  sin  ( I  ^src| /^ie) )  (11) 

where 

d,e  =  {\Rf-Rsrc?  +  zl,Y/^  (12) 

is  the  distance  from  the  source  to  the  point  (i?/,  0)  on  the  mean  plane  directly  above  the  field 
point.  As  with  the  refracted  field,  the  evanescent  field  is  found  for  a  flat  interface  making  it  a 
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zeroth-order  contribution.  If  9ie  is  greater  than  the  critical  angle,  the  evanescent  contribution  is 
set  to  zero.  Otherwise,  it  takes  the  form  [6,  see  (C25)] 


Pe(^/)  f)  \2{9ie)Psrc{i  te) / diei 


(13) 


where  Bsrc,e  is  the  source  response  function  evaluated  at  (i?/,0).  The  plane-wave  pressure 
transmission  coefficient  is 

2p  sin  9ie 


ri2(«i,)  =  i  +  r„(»„)  =  -  .  . 

P  sm  9ie  +  K  sm  92e 

where  6*2e  is  a  complex  angle  obtained  from  Snell’s  law  with  loss  neglected;  thus, 

sin6*2e  =  (1  —  COS^  9ie)^B _ 


(14) 


(15) 


Note  that  when  6*je  is  less  than  the  critical  angle,  sin  6*2e  is  purely  imaginary.  Finally,  the  complex 
time  delay  te  is 

te  =  [die  +  Kd2e)/Ci  =  {die  +  K\Zf\  Sin6'2e)/Ci.  (16) 

Insertion  of  (16)  into  (13)  reveals  the  characteristic  exponential  decay  with  depth  into  the 
sediment  of  an  evanescent  field.  Thus,  Pe{rf,t)  is  important  only  within  a  few  wavelengths 
of  the  water-sediment  interface. 

The  contribution  to  the  penetrating  acoustic  field  due  to  scattering  from  the  rough  surface  is 
obtained  with  first-order  perturbation  theory,  which  gives 


P 


d^R 


Ag{R)h{R)p 

src(t  ts 


(17) 


\rf  -  R\\R-  Vg 

(See  (Cl 5)  in  [6]  where  —kl^cl  =  — cUq  has  been  replaced  by  the  double  time  derivative  applied 
to  the  source  function.)  The  subscript  s  denotes  a  quantity  associated  with  the  scattered  field, 
while  the  subscript  2s  indicates  the  scattered  field  is  in  the  sediment  (medium  2).  The  complex 
time  delay  is  tg  =  {di  -f  Kd2g)/ci  where  di  =  |i?  —  and  d2g  =  |r/  —  R\.  The  complex 
coefficient  Ag{R)  is  given  by 

Ag{R)  =  rQBgrcTi2{di)T2l{92) 

Rj^  R  R  Rgfe 


P 


P 


X 


-f  sin  9ig  sin  6*2^ 


d2s  di 

with  the  plane-wave  pressure  transmission  coefficient  from  water  into  the  sediment 

2psm.9i 


Tl2{0^)  = 


psin  9i  +  icsin6*2s’ 


(18) 


(19) 


IEEE  JOURNAL  OE  OCEANIC  ENGINEERING,  VOL.  XX,  NO.  XX,  MARCH  2011 


and  the  coefficient  for  transmission  from  the  sediment  into  water 

2k  sin  62 

J-  — 


(20) 


ic  sin  6*2  +  p  sin  6^15 

From  the  specified  geometry,  the  real  grazing  angles  6i  and  62  satisfy  tan6*j  =  \zsrc\/\R  —  'f'srd 
and  tan  6*2  =  \zf\/\rf  —  R\,  and  the  complex  angles  623  and  9is  are  then  determined  from  Snell’s 
law: 


sin6*2s  =  (1  —  cos^  9i/ 

(21) 

sin6*is  =  (1  —  cos^  6*2)^'^^. 

(22) 

C.  Target  scattering 

The  steady-state  interaction  of  an  acoustic  field  with  a  finite  target  is  modeled  using  the 
freefield  scattering  of  an  incident  plane  wave.  At  the  ranges  commonly  found  in  SAS  applications, 
the  receiving  transducer  is  in  the  far  field  of  the  target,  and  the  scattered  pressure  for  a  freefield 
target  can  be  written  as 

Pt9t{r)  =  Pof{r,  Vi)  exp{ikr)/r,  (23) 

where  po  is  the  incident  pressure  at  the  target,  /(r,  Vi)  is  the  scattering  amplitude  where  fi 
and  f  refer  to  the  incident  and  scattering  directions,  respectively,  r  is  the  slant  range  from  the 
target  to  the  receiver,  and  k  is  the  wavenumber  in  the  water.  As  will  be  shown  later  in  this 
section,  scattering  from  a  proud  target  can  also  be  approximately  put  into  the  form  of  (23).  The 
case  of  a  buried  target  is  more  complicated,  as  discussed  in  Sec.  III-C.  The  scattering  amplitude 
also  depends  on  frequency,  but  that  dependence  is  suppressed.  For  SAS  applications  here  we  are 
interested  only  in  monostatic  backscatter,  for  which  r  =  —fi.  Restricting  the  scattering  amplitude 
to  monostatic  scattering  allows  it  to  be  expressed  using  a  single  set  of  spherical  angles  f{6,  <p), 
where  the  coordinate  origin  is  centered  on  the  target.  For  a  finite  cylindrical  target,  we  choose  9 
to  be  the  polar  angle  relative  to  the  cylinder  axis,  and  then  by  symmetry  the  scattering  amplitude 
has  no  dependence  on  ip.  For  a  spherical  target,  symmetry  demands  that  there  be  no  dependence 
on  9  or  ip.  In  what  follows  the  scattering  amplitude  is  written  as  f{9),  which  includes  both  cases. 

It  is  convenient  to  express  the  scattering  amplitude  in  terms  of  the  dimensionless  form  function 
F{9y. 

fio)  =  lF{9y 


(24) 
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where  a  is  the  radius  of  either  a  sphere  or  cylinder.  The  scattering  cross  section  per  unit  solid 
angle  is  given  by 

^  =  \fm^  =  j\m\^  =  '^\m\^-  (25) 

Thus,  the  form  function  can  be  viewed  as  an  angular  and  frequency  dependent  correction  to  the 
scattering  pattern  of  an  ideal  scatterer  with  total  cross  section  tto^  that  scatters  uniformly  over 
471  steradians.  The  target  strength  is  given  by  TS*  =  101ogio(‘^)- 

Scattering  for  spheres  and  spherical  shells  is  well  understood  and  expressions  for  various  form 
functions  are  readily  available  [8]-[ll].  In  particular,  for  a  fluid-filled,  elastic,  spherical  shell 
see  Eq.  (2)  in  [10]. 

Stanton  considered  broadside  and  near  broadside  scattering  of  sound  from  finite  fluid  and 
elastic  cylinders  [12],  [13].  In  the  simulations  of  Sec.  Ill,  the  axis  of  symmetry  of  the  finite 
cylinder  is  parallel  to  the  cross  range  axis.  Stanton’s  model  was  developed  under  the  assumptions 
that  end  effects  can  be  ignored  and  the  length  of  the  cylinder,  L,  satisfies  L  2\/r\  and  r  ^  L, 
where  A  is  the  wavelength  of  the  incident  acoustic  field  in  the  host  medium.  Stanton’s  result 
can  be  cast  into  the  following  expression  for  the  form  function  for  backscattering  from  a  finite 
cylinder: 

where  F^\k^a)  is  the  form  function  for  backscattering  of  a  normally  incident  plane  wave  from 
an  infinite  cylinder.  In  (26)  the  perpendicular  and  parallel  components  of  the  wave  vector  are 
k_i_  =  k  sin  6  and  /cy  =  k  cos  9,  respectively.  The  polar  angle  9  is  determined  by  cos  9  =  fc  -  Vgrc, 
where  fc  is  a  unit  vector  parallel  to  the  cylinder  axis  and  fgrc  is  a  unit  vector  directed  from  the 
target  center  to  the  source  position. 

To  obtain  (26)  first  recall  that  scattering  of  a  plane  wave  from  an  infinite  cylinder  is  a  two- 
dimensional  problem,  and  therefore  in  place  of  (23)  the  scattered  field  is  given  by 

Ptgt{r)  =  Pof{r,  fi)  exp{tkr)/r^^^.  (27) 

For  backscattering  from  an  infinite  circular  cylinder,  the  scattering  amplitude  has  no  angular 
dependence,  and  here  is  represented  in  terms  of  its  frequency  dependence  as  f^''{k±a).  (When 
expressed  in  terms  of  wavenumbers,  the  full  dependence  includes  the  compressional  and  shear 
wavenumbers  in  the  elastic  cylinder,  but  that  dependence  is  suppressed.)  An  expression  for 
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follows  from  Eq.  (13)  in  [13]  when  specialized  to  backscattering.  By  convention  [14] 
the  dimensionless  form  function  for  scattering  in  a  2-D  geometry  is  related  to  the  scattering 
amplitude  by  F^’'{k±a)  =  {2/ay^‘^ f^\k±a).  By  using  this  relation  and  (24)  plus  Stanton’s 
expressions  for  the  scattering  amplitudes  for  both  finite  and  infinite  length  cylinders,  (26)  readily 
follows. 

Ye  revisited  Stanton’s  model  and  provided  a  scattering  amplitude  that  includes  the  scattering 
from  a  fiat  end  [15].  His  model  was  developed  from  a  Kirchhoff  scattering  integral,  which  in  the 
backscattering  direction,  reduces  to  (26)  along  the  body  of  the  finite  cylinder.  The  backscattering 
contribution  from  a  flat  end  is  found  to  be 

Fend{d-,ka)  =  —iTii{'K/2  —  9)Ji{2kasind)  cot  6* 

X  exp(— cos  6^),  (28) 

where  an  explicit  dependence  on  ka  is  shown  to  aid  discussion  in  Sec.  III-B.  Tn  is  the  reflection 
coefficient  for  the  interface  between  the  water  and  the  flat  end,  Ji  is  the  first-order  cylindrical 
Bessel  function,  and  9  is  the  angle  that  the  incident  wave  vector  makes  with  the  cylinder  axis. 
(The  exponent  sign  and  the  overall  sign  given  in  [15]  have  been  corrected  in  (28).)  Because 
our  present  use  of  F^ndiO,  ka)  is  for  incidence  near  end-on,  we  use  for  simplicity  a  reflection 
coefficient  in  the  form  of  (3)  but  with  the  material  parameters  of  the  cylinder  and  the  surrounding 
water;  this  neglects  effects  of  shear  in  the  reflection  coefficient. 

Scattering  from  proud  targets  is  more  difficult  to  model  due  to  the  interaction  of  the  sound 
field  with  the  target  and  bottom,  including  the  possibility  of  multiple  scattering.  For  this  case 
we  have  used  a  simplified  treatment  in  which  the  interaction  of  the  incident  sound  field  with 
the  water- sediment  interface  and  the  target  is  modeled  from  a  simple  geometric  ray  analysis  [5, 
See  Fig.  5].  The  ray  diagrams  for  a  sphere  and  cylinder  at  broadside  incidence  are  identical,  and 
yield  an  approximate  backscattering  form  function 

F{9i)  =  Ttgt{Tr/2)exp{-t2ka) 

+  2rtgt{7r/2  -  9i)rii{9i)exp[t2ka{sm9i  -  cos9i)] 

+  Ttgt{7i/2)Tly9i)exp[i2ka{2sm9i  -  1)],  (29) 

where  9i  is  the  incident  grazing  angle  relative  to  the  sediment  surface.  The  form  function  given 
by  (29)  can  be  obtained  using  the  Kirchhoff  approximation  for  scattering  followed  by  a  stationary 
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phase  analysis,  consistent  with  a  geometric  ray  analysis.  (See  Baik  and  Marston  [16]  and  their 
Appendices  I  and  II;  these  authors  treat  an  analogous  problem  of  a  rigid  cylinder  near  an  interface 
and  is  a  more  general  treatment  than  required  here  in  that  the  cylinder  can  intersect  the  surface.) 
The  water-sediment  reflection  coefficient,  rii(6*i),  is  given  by  (3),  and  Ttgt  is  also  evaluated 
from  (3)  using  material  properties  of  the  target  for  the  second  medium.  The  first  term  in  (29)  is 
the  field  specularly  backscattered  from  the  target.  The  second  term  represents  a  ray  that  reflects 
from  the  water-sediment  interface  and  is  then  specularly  scattered  by  the  target  back  to  the 
source/receiver.  This  term  also  includes  the  reciprocal  path  of  first  specularly  scattering  from  the 
target  and  then  reflecting  from  the  water- sediment  interface.  The  last  term  accounts  for  a  surface 
reflected  ray  that  specularly  backscatters  from  the  target  and  then  follows  its  incoming  path  back 
to  the  source/receiver.  With  direct  substitution  of  (29)  into  (24),  the  approximate  scattering  from 
a  proud  spherical  target  is  given  by  (23).  Substitution  of  F{9i)  from  (29)  for  F^^{k±a)  in  (26) 
yields  an  approximation  for  scattering  from  a  proud  cylindrical  target  at  and  near  a  broadside 
orientation.  In  this  expression  9i  is  again  the  incident  grazing  angle  relative  to  the  sediment 
surface,  while  9  is  the  angle  the  incident  wave  vector  makes  with  the  cylinder  axis. 

Equation  (29)  includes  only  geometrically  reflected  ray  paths  and  neglects  possible  surface 
guided  waves  contributions  to  the  scattered  field  (e.g.,  leaky  Rayleigh  [5]  or  Lamb  waves  [11]). 
It  is  possible,  however,  to  include  the  most  important  elastic  effects  by  inserting  more  general 
form  functions  into  each  of  the  terms  in  (29).  The  first  step  is  to  re-express  each  term  of  (29)  in 
terms  of  the  form  function  for  scattering  from  a  sphere  or  cylinder  based  on  geometric  reflection 
only.  To  do  this  we  must  relax  our  earlier  restriction  to  only  monostatic  scattering,  since  the 
second  term  in  (29)  corresponds  to  a  vertical  bistatic  scattering  combined  with  a  reflection  from 
the  water-sediment  interface.  This  extension  of  the  formalism  will  at  first  be  restricted  here  to 
the  case  of  a  spherical  target.  Thus,  we  write  the  form  function  for  a  spherical  target  as  Faph{4>), 
where  0  is  the  vertical  angle  for  bistatic  backscatter,  defined  such  that  0  =  tt  for  monostatic 
backscatter.  In  particular,  for  scattering  based  on  geometric  reflection  (i.e.,  reflection  in  the  high 
frequency  limit)  one  obtains  Fsph{(j))  =  Ttgt{(j)/2)  exp{—i2kacos[{TT  —  4>)/2]}.  The  phase  factor 
accounts  for  the  difference  in  path  length  between  a  path  for  scattering  from  a  point  at  the  center 
of  the  sphere  and  the  path  that  reflects  from  the  surface  of  the  sphere.  [An  equivalent  phase 
factor  is  given  in  [17]  for  the  case  of  a  rigid  sphere,  Eq.  (10.103).  Note  that  (29)  corresponds 
to  using  just  the  leading  term  in  Eq.  (10.103)  for  the  spherical  target  form  function.]  Eor  the 
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geometry  of  a  proud  sphere,  0  =  tt  in  the  first  and  third  terms  in  (29)  and  0  =  tt  —  26**  in  the 
seeond  term.  Thus  (29)  becomes 

FplOi)  =  Fsphin)  +  2Tii{ei)Fsph{Ti  -  26*^)  exp(*2A;asin6'i) 

+  r^i(6'i)Fsph(7r)exp(*4A;asin6'i).  (30) 

The  phase  factors  that  survive  in  (30)  arise  from  the  path  length  differences  for  the  terms  in 
(29)  when  the  target  is  considered  a  point.  The  form  factors  in  (30)  can  now  be  generalized 
to  approximately  account  for  the  most  importance  surface  guided  wave  contributions.  This  was 
done  by  summing  the  partial  wave  series  for  scattering  from  an  elastic  spherical  shell  as  given 
by  Eq.  (3)  in  [18],  where  the  series  was  truncated  using  the  second  condition  given  in  Eq.  (16) 
in  [18]  for  all  ka.  The  partial  wave  series  gives  a  more  complete  description  of  the  specular 
reflection  component  than  given  by  (29),  and  also  includes  the  contribution  of  surface  guided 
waves.  The  latter  contribution  may  be  overestimated  because  the  effects  of  the  sediment-target 
contact  are  not  being  taken  into  account. 

The  final  piece  of  the  target  simulation  is  a  Eourier  synthesis  to  generate  a  scattered  acoustic 
pulse  emanating  from  the  target  location;  this  is  necessary  to  account  for  the  frequency  depen¬ 
dence  of  the  target’s  form  function.  Zhang  et  al.  [19]  provide  a  detailed  description  of  a  procedure 
similar  to  one  outlined  here.  The  complex  pressure,  after  propagation  from  the  source  to  the 
target  location  in  the  absence  of  the  target,  is  convolved  with  the  relevant  form  function.  This  is  a 
three  step  process.  Eirst,  the  complex  pressure  is  Eourier  transformed  into  the  frequency  domain. 
Second,  the  form  function  as  a  function  of  angular  frequency,  F{9,uj),  is  computed  for  each  of 
the  positive  frequencies  of  the  complex  signal,  and  the  negative  frequencies  are  populated  by 
F{9;  —uj)  =  F*{9;u)  where  *  denotes  complex  conjugation.  Third,  the  scattered  signal  is  given 
by  an  inverse  Eourier  transform  of  the  product  of  the  form  function  and  the  Eourier  transform 
of  the  complex  signal. 

III.  Numerical  Simulations  and  Comparisons  to  Experiments 

Several  numerical  simulations  were  conducted  to  validate  individual  models  described  in 
Sec.  II,  and  some  of  those  results  are  reported  here.  The  simplest  simulations  are  the  freefield 
scattering  from  targets.  These  permit  validation  of  the  target  models  and  the  time  domain  imaging 
(TDI)  algorithm  described  in  Appendix  II.  Proud  targets  are  then  considered,  which  introduces 
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the  interaction  of  the  SAS  signals  with  both  the  water-sediment  interface  and  target.  These 
simulations  are  compared  to  SAX04  data  and  data  collected  in  the  Naval  Surface  Warfare 
Center,  Panama  City  Division  (NSWC  PCD)  test  pond  facility.  The  simulations  culminate  with 
the  imaging  of  buried  targets  under  various  conditions  that  correspond  to  SAX04  measurements. 

Unless  otherwise  specified,  the  simulated  data  are  sampled  at  400  kHz  and  the  time  ranges 
from  0  to  24  ms.  The  length  of  the  SAS  aperture  is  typically  10  m  and  400  individual  signals 
are  generated.  This  gives  a  spatial  sampling  of  Ay  =  0.025  m,  which  is  consistent  with  the  Ay 
used  during  the  experiments. 

A.  Freefield  scattering 

For  the  freefield  simulations,  a  source  and  a  receiver  are  co-located,  and  initially  placed  at 
f'src  =  (0,  —5,  0)  m.  Omnidirectional  beam  patterns  are  used  (i.e.,  Bgrc  =  1  and  B^cv  =  !)•  A 
rigid  sphere  or  rigid  finite  cylinder  with  broadside  orientation  is  placed  at  rtgt  =  (10,  0,  0)  m 
where  the  radii  are  a  =  0.21  m  and  the  length  of  the  cylinder  is  L  =  1.1  m.  The  target  is  in 
lossless  water  with  ci  =  1482  m/s  and  pi  =  1000  kg/m^,  which  correspond  to  the  water  at  the 
NSWC  PCD  test  pond  facility  [20],  [21]. 

The  form  of  the  chirp  signal  emitted  by  the  source  is  given  in  (1.2).  The  carrier  frequency  is 
/o  =  20  kHz  with  a  bandwidth  of  fb  =  15  kHz  and  chirp  duration  of  fc  =  3  ms.  The  envelope 
parameter,  Oc  =  8,  produces  an  effective  bandwidth  of  fbe  ~  4.5  kHz  as  determined  from  the  full 
width  at  half  maximum  of  the  Fourier  spectrum  of  the  source  replica.  Throughout  the  simulations, 
a  unit  amplitude  is  assumed.  A  set  of  SAS  signals  are  now  generated  by  time  shifting  the  replica 
to  account  for  propagation  from  rgrc  to  rtgt  and  scaling  by  the  geometric  spreading,  performing 
the  convolution  described  in  Sec.  ITC,  and  finally  time  shifting  the  scattered  signal  from  Vtgt  to 
Vrcv  and  again  including  the  geometric  spreading. 

The  target  strength  of  a  sphere,  as  given  by  Urick  [22],  is  TS's  =  10  log]^o(®^/4)>  which 
yields  TSg  =  —19.58  dB.  To  determine  a  target  strength  from  Fig.  la,  it  is  noted  that  TDI 
compensates  for  the  spherical  spreading  for  propagation  from  the  source  to  each  point  in  the 
image.  The  distance  of  closest  approach  to  the  center  of  the  rigid  sphere  is  =  10  m  giving 
a  one-way  spreading  loss  of  —20.00  dB  that  needs  to  be  removed.  Additionally,  all  signals 
in  the  SAS  aperture  contribute  to  the  image  formation,  so  the  image  needs  to  be  adjusted  by 
— 201og]^Q400  =  —52.04  dB.  The  “hottest”  pixel  in  Fig.  la  corresponds  to  12.63  dB.  Thus,  we 
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Fig.  1.  SAS  images  of  the  freefield  scattering  from  (a)  a  rigid  sphere  and  (b)  a  finite  rigid  cylinder.  The  color  scale  is  referenced 
to  12.63  dB  and  11.42  dB  for  (a)  and  (b),  respectively. 


find  TSl  =  12.63  —  52.04  +  20.00  =  —19.41  dB,  which  compares  favorably  with  Urick’s  value. 

From  Uriek,  the  target  strength  of  a  finite  eylinder  at  broadside,  is  TS'c  =  10  logio(af^^/2A), 
and  yields  TSc  =  2.34  dB  at  20  kHz.  To  determine  a  target  strength  from  Fig.  lb,  the  one-way 
spreading  loss  of  —20.00  dB  again  needs  to  be  removed  to  aeeount  for  the  propagation  from  Vtgt 
to  Trcv  The  number  of  SAS  signals  that  aetively  eontribute  to  an  image  is  not  governed  by  the 
length  of  the  SAS  aperture,  but  rather  by  the  full  width  of  the  main  lobe  of  the  finite  eylinder’s 
direetivity  pattern.  From  Kinsler  et  al.  [23],  the  main  lobe  for  a  reetangular  aperture  goes  as 
sin^u/u^  with  v  =  kL  sin  9/2  and  9  is  the  angle  between  the  axis  that  is  perpendieular  to  the 
eylinder’s  axis  of  symmetry  and  the  ray  from  the  origin  to  the  field  point.  For  the  simulation 
depleted  in  Fig.  lb,  the  half-width  angle  is  9h  =  1.7°.  The  number  of  signals  eontained  within 
the  half  width  of  the  main  lobe  is  given  by  tan6*/j  =  nAy/R^  ~  9h  with  again  a  distanee  of 
elosest  approaeh  denoted  by  Rd-  The  total  number  of  eontributing  signals  is  then  N  =  2n  ^ 
29hRd/ Ay  =  23.74.  To  estimate  the  target  strength  from  Fig.  lb,  the  averaged  squared  pixel 
value  over  the  “hottest”  portion  of  the  image  is  found,  and  it  gives  101og]^Q(P^)  =  9.08  dB.  The 
estimated  target  strength  is  now  TSl  =  9-08  ~  27.51  -f  20.00  =  1.57,  whieh  is  again  in  good 
agreement  with  Uriek.  It  is  noted  that  a  value  of  77  =  21.7  reproduees  TSc. 

The  final  freefield  simulation  is  a  eement  drainage  pipe  with  an  end-on  orientation.  The 
relevant  water  parameters  from  SAX04  are  ci  =  1538  m/s  and  pi  =  1000  kg/m^.  For  the 
pipe,  the  eompressional  sound  speed  and  density  of  the  eement  are  taken  to  be  C3  =  3100  m/s 
and  p3  =  2600  kg/m^,  and  the  shear  properties  are  ignored.  The  outer  radius  of  the  pipe  is 
a  =  0.21  m  and  the  ratio  of  the  inner-to-outer  radii  is  b/a  =  0.819.  To  faeilitate  eomparison 
to  the  SAX04  measurement,  the  initial  souree  and  reeeiver  loeations  are  Vgrc  =  (0,  —5,  3.86) 
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Fig.  2.  SAS  images  of  the  freefield  scattering  from  (a)  a  cement  drainage  pipe  with  an  end-on  orientation  and  (b)  a  proud 
cement  pipe  during  SAX04. 


m  and  Vrcv  =  (0,  —5,  3.91)  m  such  that  Zsrc  and  Zrcv  are  the  nominal  measured  values  from 
the  mean  sediment  surfaee  to  the  eenter  of  the  aperture.  The  eenter  of  the  pipe  is  loeated  at 
Vtgt  =  (10.25,  0,  0.21)  m.  The  simulation  eonstruets  the  required  seattering  form  funetion  from 
two  applieations  of  (28),  F  =  Fe„d(6',  fcia)  —  F^ndid ,  kib) .  The  bandwidth  and  duration  of  the 
ehirp  have  been  inereased  to  fb  =  25  kHz  and  fc  =  4  ms  while  Oc  =  8  is  retained.  This 
gives  an  effeetive  bandwidth  of  fte  ~  7.5  kHz,  whieh  is  eonsistent  with  the  bandwidth  of  the 
SAX04  12-28  kHz  ehirp.  The  simulation  is  depieted  in  Fig.  2a  where  the  result  elearly  shows 
the  two  glints  observed  in  the  SAX04  data  shown  in  Fig.  2b.  [The  geometry  for  the  SAX04 
measurement  aetually  had  the  target  proud  on  the  sediment,  but  the  important  observation  is  that 
the  eharaeteristie  signature  of  the  target  shown  in  the  simulation  is  evident  in  the  measurement.] 
Additionally,  the  separation  distanee  between  glints  are  eonsistent  with  the  observed  glints.  Figure 
2b  also  shows  glints  at  the  far  end  of  pipe.  The  simulation  does  not  eurrently  take  into  aeeount 
the  propagation  along  the  pipe  nor  the  seattering  from  the  far  end  in  the  end-on  orientation. 

B.  Proud  target 

These  simulations  use  the  seattering  amplitude  given  by  (29)  and  (30).  A  few  eomments  are 
in  order.  First,  (29)  only  ineludes  geometrieally  seattered  eontributions  to  the  pressure,  whieh 
obviously  exeludes  possible  resonant  or  elastie  target  responses.  These  eontributions  are  ineluded 
in  (30).  Seeond,  the  simulations  also  exelude  a  surfaee  reverberation  eomponent  beeause  the 
observed  SNR  is  greater  than  25  dB.  Third,  a  set  of  signals  for  SAS  proeessing  is  generated  in 
a  manner  similar  to  the  method  deseribed  in  See.  III-A.  That  is,  a  ehirp  is  time  shifted  to  the 
target  loeation  and  sealed  by  the  propagation  distanee,  eonvolved  with  the  seattering  amplitude. 
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and  then  time  shifted  to  the  receiver.  The  following  describes  simulations  and  experiments  for 
scattering  from  a  proud  spherical  target  in  the  NSWC  PCD  test  pond  facility  and  a  cement 
drainage  pipe  at  broadside  incidence  during  SAX04. 

Figure  3  depicts  SAS  images,  generated  from  simulated  and  experimental  data,  for  the  scatter¬ 
ing  from  an  air-filled,  proud,  elastic  spherical  shell.  The  shell  has  an  outer  radius  of  a  =  0.3  m  and 
an  inner-to-outer  radii  ratio  oih/a  =  0.95.  The  shell’s  material  was  a  mild  steel.  The  SAS  aperture 
for  the  measurement  includes  764  pings,  which  gives  a  19.1  m  length  aperture.  The  center  of 
the  physical  apertures  in  the  experiment  were  initially  located  at  Vgrc  ~  (0,  —10.25,  3.86)  m  and 
Trcv  ~  (0,  —9.76,  3.89)  m.  In  the  simulation,  the  initial  source,  initial  receiver,  and  target  locations 
are  =  (0,— 5,3.86)  m,  r^cv  =  (0,— 5,3.89)  m,  and  Vtgt  =  (9.35,0.25,0.3)  m,  respectively. 
The  apertures  are  directed  toward  the  sediment  at  a  20°  depression  angle  with  source  and  receiver 
dimensions  of  0.1  x  0.05  and  0.1  x  0.1  m^.  The  chirp,  used  in  the  simulations,  has  a  carrier 
frequency  of  /o  =  20  kHz  with  a  bandwidth  of  fb  =  28  kHz,  duration  of  fc  =  4  ms,  and  envelope 
parameter  of  Oc  =  8.  This  yields  an  effective  bandwidth  of  fbe  ~  15.4  kHz,  and  the  sample  rate 
is  1  MHz.  The  density  and  sound  speed  of  the  test  pond  water  are  given  in  Sec.  III-A.  Finally,  for 
completeness,  the  density,  compressional  speed,  and  shear  speed  of  the  mild  steel  are  ps  =  7700 
kg/m^,  Qs  =  5950  m/s,  and  Cts  =  3240  m/s. 

Figures  3a  and  b  are  SAS  images  created  from  simulated  data  sets  with  the  form  functions 
given  in  (29)  and  (30).  In  (30),  we  used  the  form  function  for  an  air-filled,  elastic  spherical  shell 
for  Fsph{'K)  and  Fsph{T^  —  26*j).  Figure  3c  is  the  SAS  image  from  experimental  data.  The  three 
leftmost  features  in  these  images  can  be  uniquely  identified  with  the  geometrically  reflected  ray 
contributions  contained  in  (29).  In  Fig.  3a  and  b,  the  spatial  separation  between  the  left  and 
center  features  is  0.13  m  and  the  spatial  separation  between  the  center  and  right  features  is 
0.10  m.  For  the  experimental  data,  these  features  have  spatial  separations  of  0.085  and  0.095  m, 
respectively.  In  the  simulations,  the  target  is  assumed  to  be  in  contact  with  the  water- sediment 
interface  at  a  tangent  point.  In  the  experiment,  the  target’s  weight  slightly  depresses  the  interface 
and,  hence,  shortens  the  ray  paths  that  interact  with  the  interface.  Inspection  of  Fig.  3b  and  c 
also  reveals  a  feature  near  a  range  of  10  m.  This  feature  is  associated  with  the  subsonic  branch 
of  the  lowest  order  antisymmetric  leaky  Lamb  wave,  which  is  often  designated  as  the  Aq_  leaky 
Lamb  wave  [19]. 

The  proud  target  simulation  in  Fig.  4  approximates  a  broadside  measurement  from  a  cement 
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Fig.  3.  SAS  images  of  a  proud  elastic  spherical  shell,  (a)  The  scattering  model  used  (29)  to  generate  data.  This  simulation 
includes  only  the  geometrically  reflected  acoustic  field,  (h)  The  scattering  model  used  (30),  which  includes  both  the  reflected 
acoustic  field  and  the  elastic  response  of  the  shell,  (c)  SAS  image  generated  from  a  measurement  at  the  NSWC  PCD  test  pond 
facility. 


drainage  pipe  during  SAX04.  The  environmental  parameters  are  those  given  above  for  the  end- 
on  orientation  of  the  pipe  with  the  additional  SAX04  sediment  properties:  p2  =  2000  kg/m^, 
C2  =  1779  m/s,  and  62  =  0.01.  The  eenter  of  the  pipe  is  now  at  Vtgt  =  (9.8,0,0.21)  m,  and 
its  length  is  nominally  L  =  1.09  m.  The  loeations  and  physieal  properties  of  the  souree  and 
reeeiver  are  the  same  as  those  used  in  the  proud  spherieal  shell  measurement.  For  this  simulation 
the  form  funetion  will  again  be  generalized  to  inelude  the  most  important  elastie  effeets  using 
an  expression  similar  to  (30).  This  requires  ineluding  two  angles  in  the  finite  eylinder  form 
funetion  now  denoted  Fcy;(6*,0),  where  9  is  the  same  polar  angle  used  in  (26),  and  0  is  the 
vertieal  bistatie  angle.  The  form  funetion  given  by  (26)  is  for  monostatie  baekseatter.  In  the 
generalization  to  inelude  vertieal  bistatie  seattering  and  elastie  effeets  the  only  ehange  required 
is  to  replace  the  final  factor  by  F^’'{k±a,  (p),  corresponding  to  vertical  bistatie  scattering  from 
a  water-filled,  elastic,  infinite  cylindrical  shell  [12].  The  simulation  uses  (30)  with  the  spherical 
shell  form  factors  replaced  by  Fcyi{9,7i)  and  Fcyi{9,7i  —  29i).  At  the  extreme  ends  of  the  SAS 
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Fig.  4.  SAS  images  of  the  scattering  from  a  proud  cement  drainage  pipe  with  a  broadside  orientation:  (a)  simulation  and  (b) 
SAX04  measurement. 


aperture,  the  incident  angle  is  nominally  26°  from  broadside  and  one  may  anticipate  a  failure 
in  (30)  due  to  (26).  However,  the  discussion  of  the  freefield  image  of  the  finite  rigid  cylinder 
shown  in  Fig.  lb  suggests  that  only  the  central  ^  20  signals  contribute  significantly  to  the  image. 
The  extreme  angle  of  incidence  for  the  signals  of  consequence  is  ^  1.5°,  and  hence,  (30)  with 
(26)  should  be  adequate.  Inspection  of  Fig.  4  shows  that  the  simulation  captures  not  only  the 
length  and  width  of  the  target,  but  also  a  “doublet”.  This  splitting  is  a  clear  indication  of  the 
interaction  of  the  target  scattering  with  the  sediment  interface  as  characterized  by  the  ray  model 
in  Sec.  II-C.  Finally,  the  broad  light  blue  feature  at  an  approximate  10.2  m  range  suggests  that 
an  elastic  response  of  the  pipe  has  been  stimulated  (such  as  a  surface  elastic  wave). 

C.  Buried  targets 

The  buried  target  simulations  require  three  basic  steps:  (1)  generate  a  realization  of  the  rough 
seafloor;  (2)  compute  reverberation  from  the  seafloor;  and,  (3)  compute  a  signal  scattered  from 
the  target.  Each  step  is  briefly  outlined  below. 
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A  realization  of  a  seafloor  must  be  generated  in  accordance  with  measured  small-scale  rough¬ 
ness  and  ripple  structure.  The  small-scale  roughness  spectrum  is  given  as  an  isotropic  power  law 
with  an  algebraic  cut-off 

w{K,ky)  =  w,i{k^,  +  kiy^^^  (31) 


where  k^  and  ky  are  the  Cartesian  components  of  the  spatial  wavenumber,  k‘1  =  kl  +  ky,  and  a 
constant  cut-off  typically  taken  to  be  kc  =  20.9  m“^  (i.e.,  a  0.3  m  length  scale).  Fourier  methods 
can  be  used  to  convert  W{kx-,ky)  into  a  small-scale  surface  relief.  The  general  form  of  a  sine- 
wave  ripple  is  given  by  Jackson  [7,  see  (1)  and  (2)],  and  its  importance  to  acoustic  penetration 
is  discussed  therein.  For  simulations,  the  ripple  crests  are  parallel  to  the  cross  range  direction. 

Reverberation  is  computed  from  (1)  with  a  surface  realization  that  is  “sufficiently  large”  to 
accommodate  all  positions  along  the  SAS  aperture,  the  typically  wide  horizontal  beamwidths  of 
the  source  and  receiver,  and  the  total  duration  of  a  signal  [24].  The  surface  integral  is  converted 
to  discrete  sums  over  the  mean  surface  coordinate  Rij  =  {xi,yj)  such  that 


putn)  -^y.y. 


AttcI 

i  i=0  j=0 


(32) 


where 


=  Vr 


/? .  .  I  ir-ij  _  I  D 

-*■  1 9  •  circ.  \ 


■ij 


\rcl  Kiv  =  Arev{Rij), 

Ax  <  Ci/[4(/o  -I-  /fe/2)]  and  similarly  for  Ay,  and  =  nAt  for  n  =  0, ...,  A*.  Direct  evaluation 
of  Psrcitn  —  kpev)  for  cach  discrete  surface  point  and  discrete  round-trip  delay  time  is  possible, 
but  numerically  inefficient.  To  rapidly  evaluate  (32),  psrc{t)  is  computed  once  at  a  much  finer 
resolution  5t  At  in  the  interval  0  <  t  <  tc-  Outside  of  this  interval,  psrc{t)  is  identically 
zero.  Then,  psrdt)  is  linearly  interpolated  onto  Prevdn)  by  time  shifting  psrdt)  to  with  an 
appropriate  scaling  from  the  other  factors  in  the  summand. 

The  scattered  target  signal  is  computed  in  three  steps.  First,  the  refracted,  evanescent,  and 
surface  scattered  contributions  in  going  from  the  source  to  the  buried  target  location  are  deter¬ 
mined  from  (4),  (13),  and  (17),  respectively.  The  surface  integral  in  (17)  is  again  reduced  to  a 
double  summation  over  discrete  points  on  the  surface  realization 


P2s{tn)  = 


pAxAy  . 


4:71  d  ^  ^  yj  yj 


psrd^n  -  dd: 


(33) 


with  Mx  Nx  and  My  Ny.  Only  those  points  within  a  small  area  in  the  vicinity  of  target 
contribute.  In  (33),  rj  =  \Rij  —  rf  \  is  the  distance  from  {xi,yj)  to  the  field  point.  The  method 
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used  to  rapidly  evaluate  (32)  in  general  cannot  be  applied  to  (33)  because  fj  is  a  complex 
time  delay.  However,  for  signals  of  the  form  given  by  (L1)-(L3),  interpolation  schemes  can  be 
implement  to  account  for  fj. 

The  second  step  convolves  the  predicted  penetrating  signals  with  the  target  response  as 
described  in  Sec.  ITC.  To  investigate  the  importance  of  the  refracted,  evanescent,  and  surface 
scattered  contributions,  each  is  convolved  with  the  target  response  independently. 

The  third  step  propagates  the  three  scattered  signals  from  the  target  to  the  receiver.  To 
accomplish  this  step,  the  scattered  signals  are  again  propagated  via  (4),  (13),  and  (17)  with 
an  appropriate  change  of  coordinates  for  the  source  and  field  point  locations.  By  tracking 
the  individual  paths  back  to  the  receiver,  a  total  of  nine  possible  signals  are  generated.  For 
example,  the  acoustic  energy  refracted  into  the  sediment  scatters  from  the  target  and  a  portion 
of  this  scattered  acoustic  energy  can  arrive  at  the  receiver  via  a  refracted,  evanescent,  or  surface 
scattering  mechanism.  Hence,  the  refracted-refracted  signal  is  PRR{t),  the  refracted-evanescent 
signal  is  PRsit),  and  the  refracted-scattered  signal  is  PRi{t).  A  final  detail  of  the  third  step  is  a 
simple  time  shifting  and  interpolation  scheme  similar  to  the  one  used  in  the  surface  reverberation 
computation.  This  is  required  because  the  scattered  signals  no  longer  have  an  analytic  form.  In 
this  step,  only  attenuation  from  the  exponential  factors  in  (4),  (13),  and  (17)  is  included  while 
the  dispersion  associated  with  the  complex  time  delay  is  omitted. 

In  preparation  for  the  SAX04  SAS  measurements  conducted  during  October  2004,  proud  and 
buried  target  fields  were  deployed  during  the  Spring  of  2004  off  the  coast  of  Fort  Walton  Beach, 
Florida.  The  intent  of  the  early  deployment  was  to  allow  natural  scouring  and  benthic  activity 
to  heal  any  anthropogenic  artifacts.  Subsequent  to  the  deployment  and  prior  to  the  SAX04 
field  effort,  several  major  weather  events  nearly  destroyed  the  target  fields.  Fortunately,  a  few 
buried  targets  survived,  and  a  focusing  fluid  sphere  and  cement  drainage  pipe  were  located  from 
SAS  measurements.  For  brevity,  only  images  constructed  from  data  sequence  16  (recorded  on 
2  October  04)  and  sequences  44  and  51  (recorded  on  4  October  04)  are  analyzed  here. 

The  APL-UW  SAS  system  and  the  STMS2  electronics  package  have  been  described  elsewhere 
[4],  so  only  the  operational  parameters  relevant  to  the  SAS  processing  are  given.  An  instrumented 
tower  moved  along  the  rail  at  a  constant  speed  of  0.05  m/s,  and  it  emitted  2  pings  per  second. 
Imposing  the  “stop-start”  assumption  [25],  the  spatial  separation  of  adjacent  pings  is  Ay  =  0.025 
m.  A  typical  sequence  contains  1080  pings,  which  gives  a  27  m  SAS  aperture.  Each  ping  was 
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Fig.  5.  SAS  image  of  large  area  from  sequence  44  of  SAX04  data.  The  carrier  frequency  and  bandwidth  are  /o  =  150  and 
fb  =  80  kHz,  respectively. 


recorded  for  200  ms  at  a  1  MHz  sample  rate. 

Figure  5  displays  the  SAS  image  created  from  sequence  44.  A  110-190  kHz  chirp  with  a  4 
ms  duration  was  applied  to  the  source.  The  source  and  receiver  apertures  are  0.1  x  0.02  and 
0.1  X  0.01  m^  where  the  leading  dimension  is  the  horizontal  length.  Both  apertures  are  directed 
at  the  sediment  interface  at  a  20°  depression  angle.  The  source  and  receiver  elevations  are  4.83 
and  4.79  m  above  the  mean  sediment  interface.  Inspection  of  Fig.  5  shows  regions  of  high  and 
low  backscatter.  The  regions  of  low  backscatter  have  been  correlated  with  a  thin  layer  of  mud  on 
the  sediment  surface  via  diver  observations.  Furthermore,  measurements  of  sediment  profiles  by 
Tang  et  al.  [26]  in  nearby  locations  revealed  distinct  mud  layers.  Clearly,  if  a  buried  target  were 
located  beneath  a  mud  patch,  then  the  reduced  backscatter  would  facilitate  detection.  Although 
not  visible  due  to  the  frequency  band,  a  buried  focusing  fluid  sphere  is  detected  at  a  lower 
frequency  at  approximately  (10,10)  m,  which  is  near  the  edge  of  a  mud  patch. 

The  rough  surface  realization,  used  in  the  simulations  to  follow,  was  composed  from  two  sine- 
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wave  ripples  and  a  small-scale  isotropic  roughness  component  for  a  water-sand  interface.  The 
ripple  wavelengths  were  0.75  and  0.44  m  with  RMS  heights  of  0.015  and  0.011  m,  respectively. 
The  phases  of  the  ripples  at  a;  =  0  were  ^75  =  0°  and  044  =  37°.  These  ripples  are  consistent 
with  profiles  given  by  Tang  et  al.  [26].  The  small-scale  roughness  was  generated  from  (31) 
with  72  =  5.29  and  W2  =  1.597  which  were  obtained  by  Williams  et  al.  [27]  from 

measurements  on  a  water-sand  interface  during  SAX04.  Numerically,  the  generated  surface  was 
60  X  40  m^  with  the  first  dimension  corresponding  to  range  and  the  second  to  cross  range.  The 
spatial  increments  were  Ax  =  Ay  =  0.008  m. 

The  sound  speed  and  density  of  the  water  are  given  in  Sec.  III-A;  while  the  sediment’s  sound 
speed,  density,  and  loss  parameter  are  provided  in  Sec.  III-B.  The  sound  speeds  give  a  critical 
angle  of  30°.  The  material  properties  of  the  cement  drainage  pipe  are  given  in  Sec.  III-A.  The 
focusing  fluid  sphere  is  a  thin  steel  shell  filled  with  a  silicon  oil.  The  density,  compressional 
speed,  and  shear  speed  in  the  steel  are  ps  =  7890  kg/m^,  Q3  =  5790  m/s  and  qs  =  3100  m/s, 
respectively.  The  outer  radius  is  a  =  0.178  m  and  the  inner-to-outer  radii  ratio  is  b/a  =  0.99. 
The  properties  of  the  silicon  oil  were  selected  to  maximize  the  focusing  when  embedded  in  a 
sediment.  The  sound  speed  and  density  are  C4  =  1004  m/s  and  =  800  kg/m^.  The  silicon  oil 
is  also  given  a  small  loss  parameter  of  54  =  0.01,  which  effectively  dampens  the  elastic  response 
of  the  shell  that  is  associated  with  its  shear  properties. 

Figure  6a  displays  the  SAS  image  created  from  sequence  16.  The  source  signal  is  a  12-28 
kHz  chirp  with  a  4  ms  duration.  The  source  and  receiver  apertures  are  0.1  x  0.05  and  0.1  x  0.1 
m^,  and  the  depression  angles  are  again  at  20°.  The  source  and  receiver  elevations  were  3.86 
and  3.91  m  above  the  mean  sediment  interface,  which  corresponds  to  a  nominal  grazing  angle 
of  20°  at  the  sediment  interface  above  the  focusing  fluid  sphere.  After  the  SAS  data  collection, 
divers  probed  the  site  with  thin  rods  and  estimated  a  burial  depth  of  0.86  m  to  the  top  of  the 
shell.  Figure  6b  displays  an  image  created  from  a  simulation  of  1080  pings  based  on  the  models 
presented  here  and  our  best  estimate  of  the  environmental  conditions,  where  only  a  single  rough 
surface  realization  was  used.  The  SAS  aperture  in  the  simulation  extended  from  -13.5  to  13.5  m 
and  adjacent  pings  are  0.025  m  apart.  Each  ping  is  24  ms  in  duration.  The  center  of  the  focusing 
fluid  sphere  was  located  at  Vtgt  =  (10, 0,  —1.04)  m  (i.e.,  the  top  of  the  shell  is  buried  to  a  0.86 
m  depth  and  located  at  a  horizontal  range  of  10  m  from  the  SAS  track). 

To  achieve  a  comparable  SNR  in  Fig.  6b  to  the  observed  SNR  in  Fig.  6a,  the  reverberation 
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Fig.  6.  SAS  image  of  (a)  a  buried  focusing  fluid  sphere  from  data  sequence  16  taken  during  SAX04  (b)  a  numerical  simulation. 


predicted  by  (1)  was  reduced  by  a  multiplicative  factor  of  0.18.  At  first  glance,  a  14.9  dB 
compensation  may  suggest  that  a  problem  exists  within  the  simulation.  However,  inspection  of 
Fig.  5  indicates  that  the  overlying  mud  layer  reduced  reverberation  by  6  to  8  dB  relative  to  a 
region  without  mud.  It  is  also  noteworthy  that  roughness  parameters,  used  in  the  generation  of  the 
small-scale  surface  roughness,  are  those  for  a  sand-water  not  sand-mud  interface.  Additionally, 
the  RMS  heights  of  the  ripples  were  estimated  from  surface  profiles  taken  at  nearby  sites  not  at 
the  location  of  the  sphere.  Diver  observations  elsewhere  at  the  site  and  at  later  times  found  an 
estimated  0.07  m  peak-to-trough  ripple  amplitude.  This  suggests  that  an  increase  in  the  RMS 
ripple  heights  could  account  for  difference  in  the  observed  and  simulated  SNR.  An  increase  of 
the  RMS  heights  by  a  factor  of  2  gives  a  12  dB  increase  in  the  SNR. 

The  SAS  image  in  Fig.  7  (left)  is  created  from  sequence  51.  The  source  signal  is  a  12- 
28  kHz  chirp  with  a  4  ms  duration.  The  target  is  a  cement  drainage  pipe  in  a  nearly  end-on 
orientation.  The  source  and  receiver  elevations  were  4.72  and  4.79  m  above  the  mean  water- 
sediment  interface.  With  an  observed  horizontal  range  of  38.5  m  to  the  target,  the  incident 
grazing  angle  above  the  pipe  is  approximately  7°,  which  is  well  below  the  critical  grazing 
angle.  The  SAS  image  in  Fig.  7  (right)  is  created  from  a  simulation  of  an  finite  cylindrical 
pipe  in  an  end-on  orientation.  As  in  Sec.  III-A,  the  scattering  amplitude  is  constructed  from 
(28)  via  F  =  Fend{d,  k2a)  —  Fend{d,  k2h).  The  center  of  the  near  end  of  the  pipe  was  located  at 
'Ttgt  =  (38,0,— 0.3)  m.  The  initial  source  and  receiver  locations  are  Vgrc  =  (0,-13.5,4.72)  m 
and  Trcv  =  (0,  —13.5,4.79)  m.  Simulations  generated  1080  pings  with  Ay  =  0.025  m,  and  each 
ping  is  70  ms  in  duration. 

To  achieve  a  comparable  SNR  in  Fig.  7  to  the  observed  SNR,  the  simulated  reverberation  is 


IEEE  JOURNAL  OE  OCEANIC  ENGINEERING,  VOL.  XX,  NO.  XX,  MARCH  2011 


24 


Cross  Range  (m)  Cross  Range  (m) 

Fig.  7.  Left:  SAS  image  of  a  buried  cement  drainage  pipe  from  data  sequence  51  taken  during  SAX04.  Right:  numerical 
simulation. 


reduced  by  a  multiplicative  factor  of  0.14.  Sequence  55  is  a  data  set  for  a  110-190  kHz  chirp, 
taken  2  hours  after  sequence  51.  A  10  x  10  image,  centered  at  i?  =  (38, 10)  m  and  generated 
from  sequence  55,  reveals  that  the  observed  target  is  near  an  edge,  but  outside,  of  a  mud  patch, 
and  hence  the  17  dB  compensation  cannot  be  attributed  solely  to  reduced  surface  backscatter. 
Based  on  diver  observations,  it  was  argued  that  the  RMS  ripple  heights  could  be  increased  (to 
account  for  the  detection  in  Fig.  6a).  However,  a  17  dB  compensation  would  require  a  factor  of 
2.75  increase  (or  an  RMS  height  exceeding  0.04  m). 

IV.  Conclusion 

Models  for  reverberation,  acoustic  penetration,  and  target  scattering  have  been  combined 
into  a  unified  model  in  Sec.  II.  A  salient  feature  of  the  unified  model  is  its  incorporation 
of  the  environment  into  a  “single  realization”  modeling  capability.  The  model  can  generate 
pings  suitable  for  SAS  simulations  over  a  range  of  environmental  (and  experimental)  conditions, 
which  includes  grazing  angles  below  the  sediment’s  critical  grazing  angle.  Model  predictions 
can  then  be  compared  to  experimental  data  to  gain  insight  into  the  mechanisms  important  to 
proud  and  buried  target  detection.  Section  III  discusses  several  simulations  that  validate  the 
model  against  benchmark  problems  and  experimentally  measured  target  scattering  from  proud 
and  buried  targets. 

Although  this  article  concentrated  on  the  scattering  from  spherical  and  finite  cylindrical 
targets,  the  form  of  (23)  permits  other  target  shapes.  If  a  target’s  scattering  response  can 
be  reduced  to  a  form  function  FiO^cf)),  then  it  is  straightforward  to  investigate  the  target’s 
potential  appearance  in  SAS  images.  Methods  for  determining  F(6*,  0)  include  the  geometrical 
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theory  of  diffraction,  numerical  simulation  via  finite  elements,  and  direct  freefield  experimental 
measurement.  Obviously,  the  frequency  response  in  any  new  form  function  would  need  to  span 
the  frequency  band  of  the  source  signal. 

Two  refinements  of  the  unified  model  are  possible.  First,  the  reverberation  model  in  Sec.  IT 
A  is  restricted  to  the  backscattered  direction.  By  removing  this  restriction,  the  current  model 
would  be  applicable  to  bistatic  SAS.  Second,  at  fairly  shallow  grazing  angles  and  a  given  ripple 
wavelength,  a  portion  or  even  the  entire  frequency  band  of  a  chirp  may  exceed  an  inherent 
high-frequency  cut-off  via  the  ripple  penetration  mechanism.  To  accommodate  this  scenario,  the 
penetration  model  would  need  to  be  updated.  Two  candidates  are  the  addition  of  higher  order 
terms  from  the  perturbation  theory  that  led  to  (17)  or  possibly  a  penetration  model  based  on  the 
small- slope  approximation. 


Appendix  I 

Derivation  of  time-domain  reverberation  model 

This  appendix  provides  a  derivation  of  (1),  a  time  domain  expression  for  the  pressure  field 
backscattered  from  roughness  at  the  water-sediment  interface  where  the  sediment  is  considered  as 
a  homogeneous  absorbing  fluid.  Several  approximations  are  employed  in  obtaining  (1):  first-order 
perturbation  theory  is  used  for  the  treatment  of  scattering  from  roughness,  the  stationary  phase 
approximation  is  used  to  evaluate  the  zero-order  field  where  the  frequency  dependence  of  the 
source  and  receive  beam  patterns  is  neglected,  the  sound  speed  in  the  sediment  is  approximated 
as  independent  of  frequency,  and  the  attenuation  in  the  sediment  is  approximated  as  linearly 
increasing  with  frequency.  First-order  perturbation  theory  yields  the  scattered  field  correct  to 
first  order  in  the  small  parameter  kih,  where  ki  is  the  wave  number  in  the  water,  and  h  is  the 
rms  height  of  the  interface  roughness.  The  utility  of  making  the  final  three  approximations  is  that 
a  time  domain  result  can  be  obtained  without  the  numerical  effort  inherent  in  a  Fourier  synthesis 
approach.  The  stationary  phase  approximation  assumes  the  path  lengths  are  large  compared  to 
the  wavelength,  a  good  approximation  here.  The  assumption  that  the  sediment  sound  speed  is 
independent  of  frequency  while  the  sediment  attenuation  increases  linearly  with  frequency  entails 
a  negligible  violation  of  causality  [28].  The  derivation  will  begin  with  the  general  bistatic  case, 
and  then  specialize  to  the  case  where  the  source,  receiver,  and  scattering  patch  lie  in  the  same 
vertical  plane,  sufficient  for  the  present  applications. 
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Appendix  C  in  [6]  outlines  the  derivation  of  an  analogous  expression  for  scattering  into  the 
sediment,  the  primary  penetration  mechanism  at  low  grazing  angles,  and  the  same  approximations 
given  above  are  used  there.  In  addition,  the  relatively  narrow  bandwidths  employed  in  [6] 
allowed  the  use  of  the  baseband  representation  of  the  transmitted  pulse  with  some  simplifying 
approximations.  For  the  present  case  the  broadband  linear  chirp  signal  was  emitted  from  the 
source  at  a  center  frequency  of  20  kHz  with  a  bandwidth  of  up  to  16  kHz.  A  baseband 
representation  could  still  be  used,  but  the  simplifying  approximations  made  in  [6]  would  not 
be  appropriate.  Instead,  the  baseband  representation  has  simply  been  avoided  in  the  present 
work. 

We  let  p{r,t)  represent  the  space  and  time  dependent  pressure  field.  At  unit  distance  from 
the  source  center  on  the  transmit  beam  axis,  we  write  the  pressure  as 

Psrc(f')  Psrc^ij^^i  (f-1) 

where  psrc  is  the  maximum  rms  pressure  of  the  source  pulse  (at  pulse  center)  and  where  s{t) 
is  the  normalized  pulse  waveform  with  a  maximum  rms  value  of  unity.  (This  representation 
anticipates  the  later  stationary  phase  approximation,  and  uses  the  common  expedient  of  defining 
a  source  waveform  at  unit  distance  as  if  near  field  effects  can  be  ignored.)  Specifically, 

s{t)  =  exp{-ac[{t -ts)/ts]‘^}exp{i(j)),  (1.2) 

with 

0  =  -Uo{t  -  ts)  -  TlKc{t  -  tsf.  (1.3) 

The  pulse  is  initiated  at  t  =  0  and  in  (1.2)  the  constant  Oc  permits  some  control  over  the  width 
of  the  Gaussian  envelope  independent  of  tg.  In  (1.3)  uq  =  27r/o  with  /o  the  center  frequency, 
and  Kc  is  the  chirp  rate  defined  by  Kc  =  fb/tc  where  fb  is  the  bandwidth  and  tc  =  2ts  is  the 
duration  of  the  chirp. 

To  obtain  the  time  domain  scattered  pressure  field  we  begin  with  a  Fourier  synthesis  approach, 
and  show  later  how  that  simplifies  to  the  form  given  by  (1)  with  the  approximations  mentioned 
at  the  beginning  of  this  Appendix.  First-order  perturbation  theory  is  used  to  obtain  the  scattered 
pressure  field  in  the  frequency  domain,  which  is  denoted  by  ^|^^(r,a;),  where  the  subscript  1 
indicates  a  field  in  the  water  and  the  superscript  (1)  denotes  first-order  perturbation  theory  has 
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been  used.  The  time-domain  pressure  field  is  then  given  by 

1 

=  —  /  duj exp{—tujt), 
27r 

where 


(1.4) 


S{u)  =  /  dt  s{t)  exp{iujt). 


(1.5) 


In  order  to  obtain  the  necessary  generality,  a  self-contained  derivation  will  be  outlined  here  for 
obtaining  and  then  the  time-domain  pressure  field  using  (1.4).  The  derivation  can  be 

broken  up  into  five  steps:  (1)  Using  the  Helmholtz  integral  theorem  to  express  the  scattered  field 
in  terms  of  the  fields  just  above  and  below  the  mean  plane  of  the  water-sediment  interface  and 
in  terms  of  the  Green’s  function  for  the  case  of  a  flat  interface;  (2)  using  perturbation  theory  and 
the  boundary  conditions  to  express  the  field  quantities  at  the  mean  plane  in  terms  of  zero-order 
fields  and  h{R),  the  rough  water-sediment  interface  that  varies  with  the  horizontal  coordinate  R 
(in  this  notation,  the  3-D  position  vector  r  =  R  +  zz);  (3)  simplifying  the  resulting  expression 
for  using  integration  by  parts;  (4)  approximating  the  zero-order  field  and  the  Green’s 

function  using  the  stationary  phase  approximation;  and  (5)  performing  the  u  integration  in  (1.4) 
to  get  (1). 


A.  Getting  started  using  the  Helmholtz  integral  theorem 

To  state  the  problem,  water  above  and  sediment  below  (represented  as  a  fluid)  are  separated 
by  a  rough  interface  defined  hy  z  =  h{R)  with  the  mean  surface  at  =  0.  The  2:-axis  is  taken  as 
directed  upward  from  the  mean  plane  into  the  water,  and  the  coordinate  system  origin  is  taken 
as  in  the  mean  plane  centered  on  a  section  of  the  rough  interface  of  interest.  The  steady  state,  or 
CW,  problem  is  considered  at  this  point  for  a  particular  angular  frequency  u,  but  for  simplicity 
of  notation  the  u  dependence  will  be  suppressed  until  Sec.  TE.  A  general  source  distribution 
Ssrci'i")  is  assumed  to  reside  entirely  within  the  water  (entirely  above  the  highest  point  on  the 
rough  interface),  and  there  are  no  sources  in  the  sediment.  Because  the  source  may  have  spatial 
extent,  a  point  source  is  not  being  assumed;  in  particular,  the  source  may  have  arbitrary  beam 
patterns.  The  pressure  field  in  the  water  satisfies 


V^i/’i(r)  -f  =  - 


Ssrc{r). 


(1.6) 
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The  subscript  1  will  indicate  quantities  in  the  water  and  the  subscript  2  will  indicate  quantities 
in  the  sediment.  Thus,  ki  =  u /ci  denotes  the  wave  number  in  the  water  with  ci  the  water  sound 
speed.  The  field  that  penetrates  into  the  sediment  satisfies 

VV2(’’)  +  kl'4)2{r)  =  0,  (1.7) 


because  there  are  no  sources  there.  Also,  the  wave  number  k2  will  be  complex  because  of 
sediment  absorption.  The  fields  ipi{r)  and  satisfy  the  two-fluid  boundary  conditions  on 

=  h{R): 


^i(^)|^=MH)  =  ^2(r’)U=fc(fi),  (1-8) 

1  1 

— n  ■  Vi)i{r%=h{R)+  =  —n  ■  'V'ip2{r)\^=h{R)- ,  (1-9) 

Pi  P2 

where  pi  and  p2  are  the  densities  in  the  water  and  sediment,  respectively,  and  where  n  is  a  unit 
vector  normal  to  the  surface  h{R). 

Next  let  Vi  and  V2  be  the  volumes  in  the  regions  z  >  0  and  z  <0,  respectively,  enclosed  by 
the  two  hemispheres  that  are  formed  from  a  sphere  bisected  by  the  z  =  0  plane  and  centered 
horizontally  on  the  origin.  It  is  assumed  that  Vi  entirely  encloses  the  source  distribution  Ssrd'r’)- 
Note  that  the  hemisphere  enclosing  Vi  is  bounded  from  below  by  the  z  =  0  plane,  not  the  surface 
=  h{R).  This  choice  is  essential  for  the  development  of  perturbation  theory.  For  any  two  fields 
^(r)  and  0(r)  that  satisfy  the  Helmholtz  equation  in  Vi,  the  Helmholtz  integral  theorem  relates 
the  following  volume  integral  over  Vi  to  an  integral  over  the  surface  bounding  Vp. 


'Vr 


=  [  d^'[0(r')VV(r')  -  i/’(r')VXr')]  ■  n  , 


(1. 10) 


where  n  is  the  outwardly  directed  normal  on  bounding  Vi.  In  (1. 10)  we  now  let  i/’(r')  — > 
'^i(r')  where  we  assume  that  when  h{R)  >  0,  can  be  defined  in  the  region  0  <  2:  < 

h{R)  by  analytic  continuation.  Thus,  i/’i(r’')  is  defined  throughout  Vi.  Also  in  (1. 10)  we  let 
0(r')  — Gi{r',r),  the  Green’s  function  giving  the  field  at  r' ,  which  may  be  in  Vi  or  V2,  due 
to  a  point  source  at  r  G  Vi  in  the  presence  of  a  flat  water-sediment  interface  at  =  0.  This 
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Green’s  funetion  is  the  solution  to 


V'^Gi(r',  r)  +  kfGi^r',  r)  =  —S{r'  —  r),  r'  e  V i, 

+  klGi{r\r)  =  0,  r'  G  V2,  (Ill) 

subject  to  the  boundary  conditions 

Gi(r',r)|^,=o+  =  G'i(r',r)|y=o-  (1-12) 


and 

1  dGi{r',r)  1  dGi{r',r) 

Pi  dz'  ^,^0+  p2  dz' 

The  surface  Si  is  that  portion  of  the  ;2  =  0  plane  at  the  base  of  Vi  plus  the  hemisphere  above 
it.  By  the  usual  arguments  for  our  choices  for  ^(r)  and  0(r),  the  integral  over  the  hemispherical 
surface  — 0  as  its  radius  — >  00.  For  this  argument  to  carry  through  rigorously,  it  is  convenient  to 
assume  that  h(Ii)  goes  smoothly  to  zero  when  \R\  exceeds  some  small  fraction  of  the  hemisphere 
radius.  Even  so,  as  the  hemisphere  radius  — 00,  so  does  the  radius  of  the  region  on  the  z  =  0 


(1.13) 

z'=0- 


plane  with  roughness.  We  are  left  with  just  the  surface  integral  over  the  2: 
(1. 10)  n  =  —z.  Using  (1.6)  and  (1. 11)  in  (1. 10)  gives 


Mr)=  /  dV'Ssrc{r')Gi{r',r) 


'Vi 


dS' 


'Sr 


dz' 


'=0+ 


,  dGi{r',r) 

-  Mr) 


dz' 


z'=0+. 


When  both  r  and  r'  are  in  Vi,  the  reciprocity  condition  yields 


0  plane,  and  in 


(1.14) 


Gi{r',r)  =  Gi{r,r'), 


(1.15) 


which  is  then  used  in  (1.14)  to  replace  Gi{r',r)  with  Gi{r,r').  We  can  identify  the  first  term 
on  the  right  hand  side  in  (1.14)  as  the  zero-order  field  in  Vp. 

[  dUXc(r')Gi(r,r').  (1.16) 

Jvi 

That  is,  ^®(r)  is  the  solution  in  the  absence  of  any  roughness  at  the  water-sediment  interface 
and  is  made  up  of  a  direct  component  from  the  source  and  a  component  reflected  from  the 
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water-sediment  interface.  The  scattered  field  is  then  given  by 


-  i)f\r)  =  -  I  dS' 


'Si 


Gi{r,R 


MR 


,.dGi{r,r' 


dz' 


dz' 


'=0+J 


z'=0+ 


(1.17) 


It  will  be  seen  shortly  that  the  right  hand  side  of  (1.17)  indeed  vanishes  at  zero  order. 

To  proceed  further  the  idea  is  to  use  first-order  perturbation  theory  to  express  and  its  2:- 
derivative  on  Si  in  terms  of  the  zero-order  field  and  the  rough  surface  h{R).  The  most  convenient 
approach  to  this  end  is  to  first  develop  a  second  equation  similar  to  (1.17)  by  use  of  the  Helmholtz 
integral  theorem  in  the  volume  V2  in  the  lower  half  plane.  We  then  work  with  an  equation 
similar  to  (1. 10),  the  only  differences  being  that  ki  is  replaced  by  /c2,  the  volume  integral  is 
over  the  volume  V2  and  the  surface  integral  is  over  the  surface  S2  that  encloses  V2.  We  then 
let  >  ip2{'r’')  where  we  assume  that  when  h{R)  <  0,  can  be  defined  in  the  region 

h{R)  <  <  0  by  analytic  continuation.  Thus,  i/’2(^0  defined  throughout  V2.  We  also  let 

0(r')  — i.  Gi(r',  r)  as  before,  but  now  r'  G  V2  below  the  z  =  0  plane,  while  r  E  Vi  above  that 
plane.  As  before,  as  the  radius  of  the  hemisphere  increases  without  bound,  the  integral  over  the 
hemispherical  surface  of  5*2  0.  We  are  left  with 


0  =  /  dS' 
Jsi 


Gi{r\r)\,>=o-- 


dz' 


z’=0- 


-  M'''')\z'=o- 


dGi{r',r) 


dz' 


z'=0- 


(1.18) 


where  the  integration  is  now  over  the  entire  z  =  0  plane,  the  same  as  Si  in  (1.17). 

Terms  corresponding  to  those  on  the  left  hand  side  in  (1.17)  do  not  appear  in  (1.18):  'ip2{r) 
because  6{r'  —  r)  never  acts  in  the  integration  over  V2  and  because  there  are  no  sources 

in  V2.  There  is  a  zero-order  field  in  V2,  but  it  does  not  enter  explicitly  into  (1.18).  To  prepare 
(1.18)  for  combining  with  (1.17),  the  Green’s  functions  in  (1.18)  can  be  re-expressed  using 


Gi(r',r)|^/=o-  =  Gi(r',r)|^,=o+  =  G'i(r,  r')U'=o+ 


(1.19) 


and 


dGi{r',r) 


dz' 


p2  dGi{r',r) 


z'=0- 


pi  dz' 

P2  dGi{r,r') 


Pi 


dz' 


z'=0+ 


'=0+ 


(1.20) 
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In  (1.19)  and  (1.20)  the  first  equality  follows  from  the  boundary  eonditions  (1.12)  and  (1.13), 
and  the  seeond  equality  follows  from  the  reeiproeity  eondition  (1.15).  Using  (1.19)  and  (1.20)  in 
(1.18)  gives 


0 


Gi(r,r')U/=o+ 


Pi  diJ2{r') 
P2  dz' 
dGi{r,r') 


dz' 


h'=0- 

z'=0+- 


and  adding  (1.21)  to  (1.17)  results  in 


(1.21) 


ipiir)  - 


1 

1  d%l)2{r') 

Pi  dz' 

cx 

+ 

o 

Jl 

z'=0-- 

X  p,G,ir,R')\.  (1.22) 


From  (1.22)  it  is  finally  possible  to  see  that  the  right  hand  side  has  no  eontribution  to  zero 
order  (h{R)  =  0)),  beeause  for  that  ease  both  expressions  in  braekets  vanish  by  the  boundary 
eonditions  (1.8)  and  (1.9).  The  result  for  the  seattered  field  given  by  (1.22)  is  exaet,  exeept  for 
the  analytie  eontinuation  assumption.  Perturbation  theory  has  not  been  used. 


B.  Introducing  first-order  perturbation  theory 

The  first  step  in  bringing  in  perturbation  theory  is  to  assume  a  perturbation  expansion  of  the 
fields  above  and  below  the  z  =  0  plane: 

^i(r)  =  +  •  •  •  >  (1.23) 

'fi2{r)  =  +  •  •  •  •  (1.24) 

where  the  zero-order  fields  denoted  by  the  superseript  (0)  are  the  solutions  when  h{R)  =  0,  and 
where  the  first-order  fields  denoted  by  the  superseript  (1)  are  the  results  eorreet  to  first-order 
in  the  perturbation  theory  expansion  parameter  kih,  where  h  is  the  rms  value  of  h{R).  Thus, 
all  higher-order  terms  will  be  dropped  in  the  equations  governing  the  first-order  fields.  We  then 
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find  from  (1.22)  that  the  first-order  scattered  field  is  given  by 

v>S‘ V)  =  { [v.S‘>  (fl')  -  4"’  (fl')l 


'=0+ 


1  dip^i\r') 

1  (9i/'2^^(r') 

1 

,/=o+  P2  dz' 

- 1 

1 

O 

Jl 

X  piGi{r,B!) 


(1.25) 


To  obtain  an  expression  for  the  difference  in  the  first  set  of  brackets,  the  boundary  condition 
(1.8)  is  expanded  on  both  sides  in  a  Taylor  series  about  the  mean  plane,  and  collecting  together 
all  terms  of  first  order  yields 


=  h{R) 


- 1 

dz 

1 

o 

II 

z=0+ 

(1.26) 


Then,  the  zero-order  version  of  the  boundary  condition  (1.9)  can  be  used  to  eliminate  (r) / dz\ 
yielding 

dil)f\r) 


=  h{R)  (  ^  -  1 

\pi 


dz 


(1.27) 


z=0+ 


For  the  difference  in  the  second  set  of  brackets,  we  begin  with  the  boundary  condition  (1.9), 
use  n  =  n/\n\  with  n  =  —V±_h{R)  +  z,  and  rewrite  (1.9)  as 


Pi  L  dz 


-  VxMi^)■VxI/>l(r) 


z=h{Ry 


di)2{r) 
P2  L  dz 


-V^h{R)-V^Mr) 


(1.28) 


z=h{R)- 


where  V_l  denotes  the  horizontal  gradient  operator  (i.e.,  V  =  V_l  +  zd/dz).  Expanding  the 
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fields  in  a  Taylor’s  series  about  the  mean  plane  and  collecting  the  first-order  terms  gives 


1  c)i/'4(r) 

1  (9'i/'4(^) 

-h{R) 

Pi  dz 

z=0+ 

P2  dz 

2=0- 

X 

1 

1  d‘^ipf\r) 

Pi 

dz'^ 

isi 

II 

0 

+ 

1 

to  1 

dz"^ 

2=0- 

+  -VxMi^)■Vx^^(r)|,=o+ 

Pi 

- -V^h{R)  ■  Vxi/’f  (r)|,=o-.  (1.29) 

P2 

The  second  2:-derivatives  can  be  eliminated  using  (1.6)  and  (1.7),  and  then  all  terms  on  the  right 
hand  side  involving  can  be  re-expressed  in  terms  of  using  the  zero-order  version 

of  the  boundary  condition  (1.8).  The  result  is 


1  dilj^^\r) 
Pi  dz 


X 


1 

p2  dz 


h{R) 


<’(«)  + 

\Pi  P2J 


+ 


P2J 


(1.30) 


C.  Further  simplification  using  integration  by  parts 

When  (1.27)  and  (1.30)  are  inserted  in  (1.25),  there  are  terms  linear  in  h{R)  and  one  term 
linear  in  V±h{R).  In  this  section  integration  by  parts  is  used  to  re-express  the  term  in  V ±h{R) 
as  a  term  in  h{R),  which  is  a  useful  simplification.  Let 

^  ^j9(r)  +  4b  (^’),  (1-31) 

where  denotes  the  terms  linear  in  h{R)  and  denotes  the  term  linear  in  V ±_h{R). 

Thus,  we  have 

■<(•'0,(0  =  -  /  dS'  (7  _  1)  viA(fi')  ■ 

Jsi  \Pl  P2J 
X  piGi{r,R!). 


(1.32) 
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Integrating  (1.32)  by  parts  leads  to 


Jsi 

+  (1.33) 


No  boundary  term  appears  after  integration  by  parts  because  it  has  been  assumed  that  h{R) 
vanishes  well  before  the  boundary  is  reached.  (See  discussion  related  to  (1.14)).  The  first  term 
in  (1.33)  cancels  with  the  second  term  on  the  right  hand  side  of  (1.30)  when  it  is  inserted  into 
(1.25).  Collecting  together  the  contributions  to  (1.25)  gives 


kl  /  dS'h{R) 


'Si 


K 


l-'-]^pf\R!)Gr{r,R!) 


p  —1  dipf\r') 


kj  dz' 


dGi{r,r') 


z'=0+ 


dz' 


z'=0+ 


^V^PP{R')-V^G^{r,R') 


where  k  =  k2/ki  and  p  =  P2/ Pi- 


(1.34) 


D.  The  stationary  phase  approximation 

The  expression  given  by  (1.34)  is  valid  for  arbitrary  zero-order  fields,  and  can  be  written  in 
terms  of  wave  number  integration  over  plane  waves.  If  the  stationary  phase  approximation  is 
applied  to  such  a  representation,  the  zero-order  field  reduces  to 


=  PsrcTo 


BsTcid" )  'f'  S 


exp(*/ci|r 


.1) 


5. 


3  0rii(r’,r,^c) 


exp(*A;i|r-r;^J) 


r  —  r' 


(1.35) 


where  rg  =  1  m  and  Bsrci'f'-,  Vgrc)  is  the  source  response  function  at  field  point  r  for  the  source 
at  Tgrc-  (The  beam  pattern  is  the  square  of  the  response  function.)  Also,  is  the  location  of 
the  image  source  about  the  mean  plane,  and  rii(r,rsrc)  is  the  reflection  coefficient  for  a  ray 
from  Vsrc  to  r  via  a  reflection  from  the  mean  plane.  It  follows  from  (1.35)  that  the  zero-order 
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field  on  the  mean  plane  is  given  by 


X 


exp{zki\R'  —  Vs 


(1.36) 


|-^  ^src| 

where  Bsrc{R')  is  the  source  response  function  projected  onto  the  mean  plane,  with  the  de¬ 
pendence  of  the  source  position  suppressed,  and  where  rii(6'j)  is  the  reflection  coefficient  re¬ 
expressed  in  terms  of  the  grazing  angle  of  the  incident  ray  and  is  given  by  (3).  (The  11  subscript 
indicates  that  the  ingoing  and  outgoing  rays  are  in  medium  1.)  The  stationary  phase  approxi¬ 
mation  assumes  the  path  lengths  are  large  compared  to  the  wavelength,  a  good  approximation 
here. 

Similarly,  the  stationary  phase  approximation  applied  to  the  Green’s  function  yields 


Gl(r,i^')  =  —BrUR')[i  +  ruies)] 

4:71 


X 


exp{iki\r  —  i?'|) 


(1.37) 


\r-R\ 

where  it  is  convenient  here  to  include  the  response  function  for  the  receiver,  B^cviR'),  and  where 
9s  is  the  outgoing  grazing  angle  to  the  receiver.  The  grazing  angles  9i  and  9s  are  functions  of 
R,  though  that  dependence  is  suppressed. 

Next,  working  from  (1.35)  one  finds 


(0)/ 


dz' 


=  -tkiPsrcroBsrc{R')  sin  0* 


2'=0  + 


X  [l-rn(0,)] 


exp(^fci|i^'  -  Vsrcl) 

\R'  —  rsrrl 


(1.38) 


where  sin^j  =  Zsrc/\R'  —  'f'srcl-  Similarly, 

dGi{r,r'] 


dz' 


'=0+ 


=  -^-p-BrcviR')sm9s 
471 


X  [i-rn(0.)] 


exp(*/ci|r  —  i?'|) 


(1.39) 


\r-R\ 

In  these  expressions  the  slow  dependence  on  z'  outside  the  exponentials  was  neglected,  consistent 
with  the  stationary  phase  approximation. 
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Finally,  working  from  (1.36)  and  (1.37)  the  dot  product  of  horizontal  gradients  in  the  third 
term  of  (1.34)  becomes 

^2 

=  -^PsrcroBsrc{R')Brcv{R')  COS  6i  COS  6s  COS  0 
471 

X  [l  +  rn(0,)][l  +  rn(0s)] 


X 


exp[7A;i(  r  —  R 

'1  + 

\R'  - 

"f*  src  1 )_ 

\r-R'\\ 

R  - 

src 

(1.40) 


In  obtaining  (1.40)  the  slow  dependence  on  R  outside  the  exponentials  was  neglected.  Also,  0 
is  the  horizontal  bistatic  angle  between  the  incident  and  scattered  wave  vectors  (e.g.,  0  =  0  in 
the  forward  direction,  and  0  =  180°  for  backscatter,  or  for  vertically  bistatic  backscatter). 

Now  using  (I.36)-(I.40)  in  (1.34)  yields 

klVsTC 


= 


471 


dS'Arev{R')h{R') 


X 


exp[7/ci(  r  —  R 

:'|  + 

\R  - 

src| ) 

\r-R'\\ 

src 

(1.41) 


where  Arev{R')  is  given  by  (2)  for  bistatic  backscattering  in  the  incident  vertical  plane,  i.e., 
when  COS0  =  —1.  Equation  (1.41)  gives  the  final  CW  result,  and  it  can  now  be  used  to  obtain 
a  time  domain  result  for  the  backscattered  pressure. 


E.  The  time  domain  backscattered  pressure 


Putting  the  u  dependence  explicitly  back  into  (1.41)  gives 


0f^(r)  = 


47rcf  ^ 

X  exp{zujtr 


dS' 


Si 


Arevi^R'  ^h(^R^  ^Psrc^ 
\r  —  R'\\R'  —  Tsrrl 


(1.42) 


where  Uev  =  (k  ~  .Rl  +  ~  '>’src|)/ci  is  the  travel  time  from  the  source  to  the  point  R' 

on  the  mean  plane  and  then  to  the  field  point  at  r.  Note  that  Arev{R)  has  no  dependence 
on  u  with  the  assumptions  that  have  been  made.  It  is  evident  from  (2)  and  (3)  that  when 
the  frequency  dependence  of  the  transducer  response  functions  is  neglected,  the  only  possible 
frequency  dependence  to  A^eviR')  would  enter  through  k  =  k2/ki.  But  when  the  sound  speed  in 
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the  sediment  is  approximated  as  independent  of  frequeney,  and  the  attenuation  in  the  sediment 
is  approximated  as  linearly  inereasing  with  frequeney,  n  is  independent  of  frequeney. 

Finally,  inserting  (1.42)  into  (1.4)  and  performing  the  inverse  Fourier  transform  leads  to  (1), 
where  the  faetor  in  (1.42)  leads  to  the  double  time  derivative  of  Psrc{t)- 

Appendix  II 

Time  Domain  Imaging  Algorithm 

A  time-domain  image  is  eonstrueted  via  a  delay  and  sum  beamformer.  For  a  point  in  the 
image  plane,  loeated  at  r,  the  eoherent,  pulse-eompressed,  baseband  pressure  is 

P{'r)  =  X]  A{rsk,  Trk,  r)  exp{-tu}otk)Pk{t)  (II.l) 

k 

where  Vgk  and  r^k  are  the  souree  and  reeeiver  eoordinates  assoeiated  with  the  kth  signal. 
The  eorresponding  slant  ranges  to  the  image  point  are  =  \r  —  rsk\  and  Vrk  =  l^rk  — 
r|.  The  amplitude,  A{rsk,rrk,'r‘),  eontains  the  beam  patterns  of  the  souree  and  reeeiver  as 
well  as  eompensation  for  the  geometrie  spreading  for  the  one-way  propagation  from  Vgk  to  r. 
Contributions  from  Pk{t)  that  oeeur  outside  of  the  —6  dB  down  points  of  the  main  lobe  of  the 
beam  patterns  are  ignored.  The  exponential  faetor  in  (II.l)  eompensates  for  spherieally  diverging 
waves,  and  its  argument  ineludes  the  time  delay  to  propagate  from  the  souree  to  an  image  point 
and  then  to  the  reeeiver,  tk  =  {rgk  +  rrk)/ci.  Finally,  eo-loeation  of  the  souree  and  reeeiver  is 
not  required  by  (II.l),  so  this  beamformer  ean  be  used  for  both  monostatie  [29]  and  bistatie  SAS 
[30]  applieations. 

For  the  SAX04  SAS  measurements,  the  signals  are  real  diserete  time  signals.  In  the  simu¬ 
lations,  only  the  real  part  of  the  eomplex  signal  is  retained  beeause  a  generated  set  of  signals 
is  designed  to  simulate  a  SAX04  data  set.  Henee,  the  eomplex,  pulse-eompressed,  baseband 
pressure,  Pk{t),  is  obtained  by  applying  a  Hilbert  transform  to  pk{t),  a  matched  filter,  and 
removal  of  the  carrier  frequency.  For  SAS  processing  of  the  simulated  data,  the  matched  filter 
is  constructed  from  a  Hilbert  transformed  replica  of  a  unit- amplitude  transmitted  pulse.  For  the 
SAX04  data,  the  matched  filter  is  constructed  from  a  Hilbert  transform  of  the  signal  applied  to 
the  source  transducer.  The  match  filter  is  then  normalized  to  give  a  pulse-compressed  baseband 
replica  with  an  amplitude  of  2. 
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To  form  an  image,  the  diserete  representation  of  the  image  coordinates  are 


Xi  =  {xq  -  xi/2)  +  {i -1),  (*  =  !,•■ 

■■  .IVJ. 

(II.2) 

O 

1 

+ 

1 

••  w,). 

(II.3) 

The  center  of  the  image  is  at  {xo,yo),  and  the  lengths  of  the  sides  of  the  image  are  xi  and  yi. 
Additionally,  the  image  plane  is  not  confined  to  the  Zq  =  0  plane;  however,  zq  is  assumed  to  be 
constant.  For  the  kth  signal,  the  slant  ranges  from  the  source  and  receiver  to  {xi,yj)  are 

Tskij  =  [{Xi  -  Xskf  +  {yj  -  yskY  +  {Zq  -  ZskfY'‘^,  (11.4) 

^rkij  [(^*  Xfk^  T  (l/j  Vrk^  T  (^0  ^rfc)  ]  ^  •  (II'5) 

Inspection  of  (II.l)  and  the  discrete  nature  of  the  sampled  time  signals,  Pkitn),  (n  =  I,  -  ■  ■  ,  A)) 
suggests  that  tk  seldom  coincides  with  a  discrete  time  point  Thus,  the  summation  in  (II.l) 
requires  interpolation  of  the  signals.  The  interpolation  has  been  performed  with  2,  4,  6,  and 
8-point  Lagrange  interpolation  algorithms.  Cubic  spline  interpolation  with  6,  12,  24,  and  48 
points,  constructed  on  intervals  centered  about  a  tn,  were  also  investigated.  Although  cubic 
spline  interpolation  permitted  the  construction  of  an  image,  4-point  Lagrange  interpolation  was 
found  to  be  sufficient  for  image  formation  and  is  more  numerically  efficient  than  cubic  splines. 
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